Tables and Figures

library(ggplot2)
library(RColorBrewer)
library(magrittr)
library(wesanderson)


load("../Data/Results.RData")

Data

Wilcox Test

PA_wilcox

Haplotype Statistics

HBStat_sum

Prediction Accuracy

CV_sum

FV2 Composition

FV2Comp

Tables

Mantel correlation

Man_mat
##         FixedHB HaploView HaploBlocker
## DH_KE 0.9955470 0.9800637    0.9298979
## DH_PE 0.9956359 0.9798048    0.9411741
## GC_KE 0.9921593 0.9657868    0.9440275
## GC_PE 0.9953509 0.9717126    0.9842589

Wilcoxon test of prediction accuracy (scenario1)

sign_mat <- matrix(nrow = 3, ncol = 3,
                   dimnames = list(c("FixedHB", "HaploView", "HaploBlocker"),
                                   c("Decreased", "Increased", "Equal")))
for (m in c("FixedHB", "HaploView", "HaploBlocker")){
  
  PA_wilcox_sub <- subset(PA_wilcox, HBMethod == m )
  
  sign_mat[ m, "Decreased"] <- with(PA_wilcox_sub,
                                   sum(Difference < 0 & pvalue_bon < 0.05))
  sign_mat[ m, "Increased"] <- with(PA_wilcox_sub,
                                   sum(Difference > 0 & pvalue_bon < 0.05))
  sign_mat[ m, "Equal"] <- with(PA_wilcox_sub,
                                   sum(pvalue_bon > 0.05))
}

sign_mat
##              Decreased Increased Equal
## FixedHB              0        10    10
## HaploView            0        10    10
## HaploBlocker         4         6    10

Number of haplotypes (dimention reduction)

nHB_Table <- cbind(
  subset(HBStat_sum, select = nHB_filter,
       Scenario == "Scenario1" &
         HBMethod == "FixedHB" &
         Param1 == "win20"),
subset(HBStat_sum, select = nHB_filter,
       Scenario == "Scenario1" &
         HBMethod == "HaploView" &
         Param1 == "GAB"),
subset(HBStat_sum, select = nHB_filter,
       Scenario == "Scenario1" &
         HBMethod == "HaploBlocker" &
         Param1 == "win20" &
         Param2 == "tc99")
) %>%
  set_colnames(c("FixedHB", "HaploView", "HaploBlocker")) %>%
  set_rownames(c("DH_KE", "DH_PE", "GC_KE", "GC_PE"))

print(nHB_Table)
##       FixedHB HaploView HaploBlocker
## DH_KE   84406     51860        10225
## DH_PE  106620     77314        14236
## GC_KE  112092     72215        53178
## GC_PE  127420     93445        41239

Figures

Haplotype Charateristics

# Subset Data #
dat_plot <- subset(HBStat_sum,
                   Scenario == "Scenario1")

# Plot #
ggplot(data = subset(dat_plot, HBMethod == "FixedHB"),
         aes(x = nMR_Mean, y = nHB_filter/1000, group = Train, color = Train)) +
    geom_point(shape = 4, size = 4) +
    geom_line() +
    theme_bw(base_size = 20) +
    scale_x_continuous(breaks = c(5,10,20,50,100)) +
    labs(title = "FixedHB",
         x = "Haplotype Length (SNPs)",
         y = "Total Number of Alleles (K)",
         color = "Population") +
    scale_color_manual(values = c(wes_palette(name="Darjeeling1")[c(1,3,5)],
                                  wes_palette(name="Rushmore")[c(3)])) 

# Subset Data #
dat_plot <- subset(HBStat_sum,
                   Scenario == "Scenario1")

# Plot #
  ggplot(data = subset(dat_plot, HBMethod == "HaploView"),
         aes(x = nMR_Mean, y = nHB_filter/1000, group = Train, color = Train, shape = Param1)) +
    geom_point(size = 4) +
    labs(title = "HaploView",
         x = "Haplotype Length (SNPs)",
         y = "Total Number of Alleles (K)",
         color = "Population",
         shape = "Algorithm") +
    theme_bw(base_size = 20) +
    scale_shape_manual(values = 16:18,
                       labels = c("4GAM", "GABRIEL", "SPINE")) +
    scale_color_manual(values = c(wes_palette(name="Darjeeling1")[c(1,3,5)],
                                  wes_palette(name="Rushmore")[c(3)])) 

  # Subset Data #
dat_plot <- subset(HBStat_sum,
                   Scenario == "Scenario1")

# Plot #
ggplot(data = subset(dat_plot, HBMethod == "HaploBlocker"),
         aes(x = nMR_Mean, y = nHB_filter/1000, group = Train, color = Train)) +
    geom_point(size = 3) +
    labs(title = "HaploBlocker",
         x = "Haplotype Length (SNPs)",
         y = "Total Number of Alleles (K)",
         color = "Population") +
    theme_bw(base_size = 20) +
    scale_color_manual(values = c(wes_palette(name="Darjeeling1")[c(1,3,5)],
                                  wes_palette(name="Rushmore")[c(3)]))  

Prediction Accuracy (SNP vs HB)

# Scenario 1 #
# Select the optimal haplotype sets #
sel <- with(CV_sum, (Scenario == "Scenario1") &
              ((HBMethod == "FixedHB" & Param1 == "win20")|
                 (HBMethod == "HaploView" & Param1 == "GAB")|
                 (HBMethod == "HaploBlocker" & Param1 == "win20" & Param2 == "tc99"))
            )

dat_plot <- CV_sum[ sel,]
dat_plot_snp <- subset(CV_sum,
                       Scenario == "Scenario1" & HBMethod == "SNP")

key <- match(paste(dat_plot$Trait, dat_plot$Train),
             paste(dat_plot_snp$Trait, dat_plot_snp$Train))
dat_plot$PA_SNP <- dat_plot_snp[key, "PA"]
dat_plot$HBMethod <- factor(dat_plot$HBMethod,
                              levels = c("FixedHB", "HaploView", "HaploBlocker"))
  

# Fig_scen1_SNPvsHB #  
ggplot(data = dat_plot, aes(x = PA_SNP, y = PA, color = Valid, shape = Trait)) +
    geom_abline(slope = 1, intercept = 0) +
    geom_point(size = 3) +
    scale_shape_manual(values = c(8, 15:18)) +
    facet_grid(cols = vars(HBMethod)) +
    theme_bw(base_size = 20) +
    labs(title = "Scenario 1",
         x = "Prediction Accuracy (SNP)",
         y = "Prediction Accuracy (HB)",
         color = "Population", shape = "Trait") +
    scale_color_manual(values = c(wes_palette(name="Darjeeling1")[c(1,3,5)],
                                  wes_palette(name="Rushmore")[c(3)]))    

# Scenario 2 #
# Select the optimal haplotype sets #
sel <- with(CV_sum, (Scenario == "Scenario2") &
              ((HBMethod == "FixedHB" & Param1 == "win20")|
                 (HBMethod == "HaploView" & Param1 == "GAB")|
                 (HBMethod == "HaploBlocker" & grepl(Train, pattern = "KE") &
                    Param1 == "MinSubgroup40" & Param2 == "MCMB1")|
                 (HBMethod == "HaploBlocker" & grepl(Train, pattern = "PE") &
                    Param1 == "MinSubgroup40" & Param2 == "MCMB1250"))
            )

dat_plot <- CV_sum[ sel,]
dat_plot_snp <- subset(CV_sum, 
                       Scenario == "Scenario2" & HBMethod == "SNP")

key <- match(paste(dat_plot$Trait, dat_plot$Train),
             paste(dat_plot_snp$Trait, dat_plot_snp$Train))
dat_plot$PA_SNP <- dat_plot_snp[key, "PA"]
dat_plot$HBMethod <- factor(dat_plot$HBMethod,
                              levels = c("FixedHB", "HaploView", "HaploBlocker"))
  

# Fig_scen2_SNPvsHB #  
ggplot(data = dat_plot, aes(x = PA_SNP, y = PA, color = Valid, shape = Trait)) +
    geom_abline(slope = 1, intercept = 0) +
    geom_point(size = 3) +
    scale_shape_manual(values = c(8, 15:18)) +
    facet_grid(cols = vars(HBMethod)) +
    theme_bw(base_size = 20) +
    labs(title = "Scenario 2",
         x = "Prediction Accuracy (SNP)",
         y = "Prediction Accuracy (HB)",
         color = "Population", shape = "Trait") +
    scale_color_manual(values = c(wes_palette(name="Darjeeling1")[c(1,3,5)],
                                  wes_palette(name="Rushmore")[c(3)]))    

# Scenario 3 #
# Select the optimal haplotype sets #
sel <- with(CV_sum, (Scenario == "Scenario3") &
              ((HBMethod == "FixedHB" & Param1 == "win10")|
                 (HBMethod == "HaploView" & Param1 == "GAB")|
                 (HBMethod == "HaploBlocker" & grepl(Train, pattern = "DH") &
                    Param1 == "MinSubgroup20" & Param2 == "MCMB1250")|
                 (HBMethod == "HaploBlocker" & grepl(Train, pattern = "GC") &
                    Param1 == "MinSubgroup0" & Param2 == "MCMB1"))
            )

dat_plot <- CV_sum[ sel,]
dat_plot_snp <- subset(CV_sum, 
                       Scenario == "Scenario3" & HBMethod == "SNP")

key <- match(paste(dat_plot$Trait, dat_plot$Train),
             paste(dat_plot_snp$Trait, dat_plot_snp$Train))
dat_plot$PA_SNP <- dat_plot_snp[key, "PA"]
dat_plot$HBMethod <- factor(dat_plot$HBMethod,
                              levels = c("FixedHB", "HaploView", "HaploBlocker"))
dat_plot$DHorGC <- substr(dat_plot$Train, start = 1, stop = 2)


# Fig_scen3_SNPvsHB #  
ggplot(data = dat_plot, aes(x = PA_SNP, y = PA, color = Valid, shape = Trait)) +
      geom_abline(slope = 1, intercept = 0) +
      geom_point(size = 4) +
      scale_shape_manual(values = c(8, 15:18)) +
      facet_grid(rows = vars(DHorGC), cols = vars(HBMethod)) +
      theme_bw(base_size = 20) +
      labs(title = "Scenario 3",
           x = "Prediction Accuracy (SNP)",
           y = "Prediction Accuracy (HB)",
           color = "Prediction Set", shape = "Trait")+
      scale_color_manual(values = c(wes_palette(name="Darjeeling1")[c(1,3,5)],
                                    wes_palette(name="Rushmore")[c(3)]))    

Impact of Parameters

Scenario 1

HaploBlocker

dat_plot <- subset(HBStat_sum,
                   Scenario == "Scenario1" &
                     HBMethod == "HaploBlocker" &
                     !Param2 == "NoTC")

dat_plot$Param1 <- factor(dat_plot$Param1,
                          levels = c("win0", "win5", "win10", "win20", "win50", "win100" 
                                     ),
                          labels = c("MultiWin", "5", "10", "20", "50", "100"))
dat_plot$Param2 <- factor(dat_plot$Param2,
                          levels = c("tc80", "tc85", "tc90", "tc95", "tc99"),
                          labels = c("80", "85", "90", "95", "99"))

cols <- c("MultiWin" = "#4DAF4A", # green
              "5" = brewer.pal(8, "RdYlBu")[8],
              "10" = brewer.pal(8, "RdYlBu")[7],
              "20" = brewer.pal(8, "RdYlBu")[6],
              "50" = brewer.pal(8, "RdYlBu")[2],
              "100" = brewer.pal(8, "RdYlBu")[1]
              )

## Params vs Length (HaploBlocker) ##
ggplot(data = dat_plot, aes(x = Param2, y = nMR_Mean, color = Param1, group = paste(Param1, Train))) +
      geom_line() +
      geom_point() +
      labs(title = "HaploBlocker (scenario 1)",
           x = "Target Coverage",
           y = "Haplotype Length (SNPs)",
           color = "Window Size") +
      facet_grid(rows = vars(Train), scale = "free") +
      theme_bw(base_size = 20) +
      theme(axis.text.x=element_text(angle = -45, hjust = 0)) +
      scale_color_manual(values = cols)

# Data preparation #
dat_plot <- subset(CV_sum,
                   Scenario == "Scenario1" &
                     HBMethod == "HaploBlocker" &
                     !Param2 == "NoTC")

dat_plot <- aggregate(data = dat_plot,
          PA ~ Train + Param1 + Param2,
          mean)

dat_plot$Param1 <- factor(dat_plot$Param1,
                          levels = c("win0", "win5", "win10", "win20", "win50", "win100" 
                                     ),
                          labels = c("MultiWin", "5", "10", "20", "50", "100"))
dat_plot$Param2 <- factor(dat_plot$Param2,
                          levels = c("tc80", "tc85", "tc90", "tc95", "tc99"),
                          labels = c("80", "85", "90", "95", "99"))


## Params vs PA (HaploBlocker) ##
ggplot(data = dat_plot, aes(x = Param2, y = PA, color = Param1, group = paste(Param1, Train))) +
      geom_line() +
      geom_point() +
      labs(title = "HaploBlocker (scenario 1)",
           x = "Target Coverage",
           y = "Prediction Accuracy",
           color = "Window Size") +
      facet_grid(rows = vars(Train), scales = "free") +
      theme_bw(base_size = 20) +
      theme(axis.text.x=element_text(angle = -45, hjust = 0)) +
      scale_color_manual(values = cols)

FixedHB

# Data preparation #
dat_plot <- subset(CV_sum,
                   Scenario == "Scenario1" & HBMethod == "FixedHB")
dat_plot_snp <- subset(CV_sum,
                   Scenario == "Scenario1" & HBMethod == "SNP")
dat_plot$Param1 <- factor(dat_plot$Param1,
                          levels = c("win5", "win10", "win20", "win50", "win100"),
                          labels = c("5", "10", "20", "50", "100"))
  
key <- match(paste(dat_plot$Trait, dat_plot$Train),
             paste(dat_plot_snp$Trait, dat_plot_snp$Train))


## Params vs PA (FixedHB) ##
ggplot(data = dat_plot,
       aes(x = Param1, y = PA, color = Trait, shape = Trait, group = paste(Trait, Train))) +
    geom_point(size = 2) +
    geom_line() +
    geom_abline(data = dat_plot_snp, aes(intercept = PA, slope = 0, color = Trait), linetype = "dashed") +
    scale_shape_manual(values=c(8, 15:18)) +
    facet_grid(cols = vars(Train), scales = "free_y") +
    labs(title = "FixedHB (scenario 1)",
         x = "Window Size",
         y = "Prediction Accuracy") +
    theme_bw(base_size = 20) +
    theme(axis.text.x=element_text(angle = -45, hjust = 0)) +
    guides(color = guide_legend("Trait"), shape = guide_legend("Trait")) +
    scale_color_brewer(palette = "Set1") 

HaploView

# Data preparation #
dat_plot <- subset(CV_sum,
                   Scenario == "Scenario1" & HBMethod == "HaploView")
dat_plot_snp <- subset(CV_sum,
                   Scenario == "Scenario1" & HBMethod == "SNP")
  
dat_plot <- merge(dat_plot, dat_plot_snp,
        by = c("Trait", "Train"),
        suffixes = c("_HB", "_SNP"))
  
dat_plot$PA_diff <- dat_plot$PA_HB - dat_plot$PA_SNP  
  

## Params vs PA (HaploView) ##
ggplot(data = dat_plot,
           aes(x = Trait, y = PA_diff, fill = Param1_HB, group = Param1_HB)) +
    geom_bar(position="dodge", stat="identity") +
    facet_grid(cols = vars(Train), scales = "free_y") +
    labs(title = "HaploView (scenario 1)",
         x = "Methods",
         y = "Relative Prediction Accuracy") +
    theme_bw(base_size = 20) +
    scale_fill_brewer(palette = "Set1", name = "Methods") + 
    theme(axis.text.x=element_text(angle = -45, hjust = 0))

Scenario 2

HaploBlocker

dat_plot <- subset(HBStat_sum,
                   Scenario == "Scenario2" &
                     HBMethod == "HaploBlocker")

dat_plot$Param1 <- factor(dat_plot$Param1,
                          levels = c("MinSubgroup0", "MinSubgroup5","MinSubgroup20", "MinSubgroup40", "MinSubgroup80", "MinSubgroup160"),
                          labels = c("0", "5","20", "40", "80", "160"))

dat_plot$Param2 <- factor(dat_plot$Param2,
                          levels = c("MCMB1", "MCMB1250", "MCMB5000", "MCMB20000", "MCMB80000"),
                          labels = c("1", "1250","5000", "20000", "80000"))


ggplot(data = dat_plot,
       aes(y = nMR_Mean, x = Param2, group = paste(Train, Param1),
               shape = Param1, color = Train))+
    geom_point(size = 4) +
    geom_line() +
    labs(title = "HaploBlocker Length (scenario 2)",
         x = "MCMB",
         y = "Haplotype Length",
         color = "Landrace",
         shape = "Min Subgroup") +
    theme_bw(base_size = 20) +
    scale_color_brewer(palette = "Set1") +
    guides(color = guide_legend(order = 1)) +
    theme(axis.text.x=element_text(angle = -45, hjust = 0))

dat_plot <- subset(HBStat_sum,
                   Scenario == "Scenario2" &
                     HBMethod == "HaploBlocker")

dat_plot$Param1 <- factor(dat_plot$Param1,
                          levels = c("MinSubgroup0", "MinSubgroup5","MinSubgroup20", "MinSubgroup40", "MinSubgroup80", "MinSubgroup160"),
                          labels = c("0", "5","20", "40", "80", "160"))

dat_plot$Param2 <- factor(dat_plot$Param2,
                          levels = c("MCMB1", "MCMB1250", "MCMB5000", "MCMB20000", "MCMB80000"),
                          labels = c("1", "1250","5000", "20000", "80000"))


ggplot(data = dat_plot,
       aes(y = GenomeCoverage, x = Param2, group = paste(Train, Param1),
               shape = Param1, color = Train))+
    geom_point(size = 4) +
    geom_line() +
    labs(title = "HaploBlocker Coverage (scenario 2)",
         x = "MCMB",
         y = "Genome Coverage",
         color = "Landrace",
         shape = "Min Subgroup") +
    theme_bw(base_size = 20) +
    scale_color_brewer(palette = "Set1") +
    guides(color = guide_legend(order = 1)) +
    theme(axis.text.x=element_text(angle = -45, hjust = 0))

dat_sub <- subset(CV_sum,
                   Scenario == "Scenario2" &
                     HBMethod == "HaploBlocker" &
                    !Param1 %in% c("MinSubgroup30", "MinSubgroup50"))

dat_sub$Param1 <- factor(dat_sub$Param1,
                          levels = c("MinSubgroup0", "MinSubgroup5","MinSubgroup20", "MinSubgroup40", "MinSubgroup80", "MinSubgroup160"),
                          labels = c("0", "5","20", "40", "80", "160"))

dat_sub$Param2 <- factor(dat_sub$Param2,
                          levels = c("MCMB1", "MCMB1250", "MCMB5000", "MCMB20000", "MCMB80000"),
                          labels = c("1", "1250","5000", "20000", "80000"))


dat_sub$Pop <- substr(dat_sub$Train, start = 4, stop = 5)

CV_merge <- aggregate(data = dat_sub,
                      PA ~ Pop + Param1 + Param2,
                      mean)

dat_plot <- with(data = CV_merge,
                   data.frame(PA = PA,
                              Param = c(Param1, Param2),
                              Pop = Pop,
                              Group = rep(c("subpop", "MCMB"), each = nrow(CV_merge))))
    
dat_plot$Group <- factor(dat_plot$Group,
                           levels = c("MCMB", "subpop"),
                           labels = c("MCMB", "Min Subgroup"))


  
ggplot(data = dat_plot, aes(x = Param, y = PA)) +
    geom_boxplot() +
    facet_grid(cols = vars(Group), rows = vars(Pop), scales = c("free")) +
    labs(title = "HaploBlocker Optimization (scenario 2)",
         x = "Parameters",
         y = "Prediction Accuracy") +
    theme_bw(base_size = 20) +
    theme(axis.text.x=element_text(angle = -45, hjust = 0))  

FixedHB

# Data preparation #
dat_plot <- subset(CV_sum,
                   Scenario == "Scenario2" & HBMethod == "FixedHB")
dat_plot_snp <- subset(CV_sum,
                   Scenario == "Scenario2" & HBMethod == "SNP")
dat_plot$Param1 <- factor(dat_plot$Param1,
                          levels = c("win5", "win10", "win20", "win50", "win100"),
                          labels = c("5", "10", "20", "50", "100"))
  
key <- match(paste(dat_plot$Trait, dat_plot$Train),
             paste(dat_plot_snp$Trait, dat_plot_snp$Train))


## Params vs PA (FixedHB) ##
ggplot(data = dat_plot,
       aes(x = Param1, y = PA, color = Trait, shape = Trait, group = paste(Trait, Train))) +
    geom_point(size = 2) +
    geom_line() +
    geom_abline(data = dat_plot_snp, aes(intercept = PA, slope = 0, color = Trait), linetype = "dashed") +
    scale_shape_manual(values=c(8, 15:18)) +
    facet_grid(cols = vars(Train), scales = "free_y") +
    labs(title = "FixedHB (scenario 2)",
         x = "Window Size",
         y = "Prediction Accuracy") +
    theme_bw(base_size = 20) +
    theme(axis.text.x=element_text(angle = -45, hjust = 0)) +
    guides(color = guide_legend("Trait"), shape = guide_legend("Trait")) +
    scale_color_brewer(palette = "Set1") 

Scenario 3

HaploBlocker

dat_plot <- subset(HBStat_sum,
                   Scenario == "Scenario3" &
                     HBMethod == "HaploBlocker")

dat_plot$Param1 <- factor(dat_plot$Param1,
                          levels = c("MinSubgroup0", "MinSubgroup5","MinSubgroup20", "MinSubgroup40", "MinSubgroup80", "MinSubgroup160"),
                          labels = c("0", "5","20", "40", "80", "160"))

dat_plot$Param2 <- factor(dat_plot$Param2,
                          levels = c("MCMB1", "MCMB1250", "MCMB5000", "MCMB20000", "MCMB80000"),
                          labels = c("1", "1250","5000", "20000", "80000"))


ggplot(data = dat_plot,
       aes(y = nMR_Mean, x = Param2, group = paste(Train, Param1),
               shape = Param1, color = Train))+
    geom_point(size = 4) +
    geom_line() +
    labs(title = "HaploBlocker Length (scenario 3)",
         x = "MCMB",
         y = "Haplotype Length",
         color = "Population",
         shape = "Min Subgroup") +
    theme_bw(base_size = 20) +
    scale_color_brewer(palette = "Set1") +
    guides(color = guide_legend(order = 1)) +
    theme(axis.text.x=element_text(angle = -45, hjust = 0))

dat_plot <- subset(HBStat_sum,
                   Scenario == "Scenario3" &
                     HBMethod == "HaploBlocker")

dat_plot$Param1 <- factor(dat_plot$Param1,
                          levels = c("MinSubgroup0", "MinSubgroup5","MinSubgroup20", "MinSubgroup40", "MinSubgroup80", "MinSubgroup160"),
                          labels = c("0", "5","20", "40", "80", "160"))

dat_plot$Param2 <- factor(dat_plot$Param2,
                          levels = c("MCMB1", "MCMB1250", "MCMB5000", "MCMB20000", "MCMB80000"),
                          labels = c("1", "1250","5000", "20000", "80000"))


ggplot(data = dat_plot,
       aes(y = GenomeCoverage, x = Param2, group = paste(Train, Param1),
               shape = Param1, color = Train))+
    geom_point(size = 4) +
    geom_line() +
    labs(title = "HaploBlocker Coverage (scenario 3)",
         x = "MCMB",
         y = "Genome Coverage",
         color = "Population",
         shape = "Min Subgroup") +
    theme_bw(base_size = 20) +
    scale_color_brewer(palette = "Set1") +
    guides(color = guide_legend(order = 1)) +
    theme(axis.text.x=element_text(angle = -45, hjust = 0))

dat_sub <- subset(CV_sum,
                   Scenario == "Scenario3" &
                     HBMethod == "HaploBlocker")

dat_sub$Param1 <- factor(dat_sub$Param1,
                          levels = c("MinSubgroup0", "MinSubgroup5","MinSubgroup20", "MinSubgroup40", "MinSubgroup80", "MinSubgroup160"),
                          labels = c("0", "5","20", "40", "80", "160"))

dat_sub$Param2 <- factor(dat_sub$Param2,
                          levels = c("MCMB1", "MCMB1250", "MCMB5000", "MCMB20000", "MCMB80000"),
                          labels = c("1", "1250","5000", "20000", "80000"))


dat_sub$Pop <- substr(dat_sub$Train, start = 1, stop = 2)

CV_merge <- aggregate(data = dat_sub,
                      PA ~ Pop + Param1 + Param2,
                      mean)

dat_plot <- with(data = CV_merge,
                   data.frame(PA = PA,
                              Param = c(Param1, Param2),
                              Pop = Pop,
                              Group = rep(c("subpop", "MCMB"), each = nrow(CV_merge))))
    
dat_plot$Group <- factor(dat_plot$Group,
                           levels = c("MCMB", "subpop"),
                           labels = c("MCMB", "Min Subgroup"))


  
ggplot(data = dat_plot, aes(x = Param, y = PA)) +
    geom_boxplot() +
    facet_grid(cols = vars(Group), rows = vars(Pop), scales = c("free")) +
    labs(title = "HaploBlocker Optimization (scenario 3)",
         x = "Parameters",
         y = "Prediction Accuracy") +
    theme_bw(base_size = 20) +
    theme(axis.text.x=element_text(angle = -45, hjust = 0))  

FixedHB

# Data preparation #
dat_plot <- subset(CV_sum,
                   Scenario == "Scenario3" & HBMethod == "FixedHB")
dat_plot_snp <- subset(CV_sum,
                   Scenario == "Scenario3" & HBMethod == "SNP")
dat_plot$Param1 <- factor(dat_plot$Param1,
                          levels = c("win5", "win10", "win20", "win50", "win100"),
                          labels = c("5", "10", "20", "50", "100"))
  
key <- match(paste(dat_plot$Trait, dat_plot$Train),
             paste(dat_plot_snp$Trait, dat_plot_snp$Train))


## Params vs PA (FixedHB) ##
ggplot(data = dat_plot,
       aes(x = Param1, y = PA, color = Trait, shape = Trait, group = paste(Trait, Train))) +
    geom_point(size = 2) +
    geom_line() +
    geom_abline(data = dat_plot_snp, aes(intercept = PA, slope = 0, color = Trait), linetype = "dashed") +
    scale_shape_manual(values=c(8, 15:18)) +
    facet_grid(cols = vars(Train), scales = "free_y") +
    labs(title = "FixedHB (scenario 3)",
         x = "Window Size",
         y = "Prediction Accuracy") +
    theme_bw(base_size = 20) +
    theme(axis.text.x=element_text(angle = -45, hjust = 0)) +
    guides(color = guide_legend("Trait"), shape = guide_legend("Trait")) +
    scale_color_brewer(palette = "Set1") 

FV2 Component

par(mfrow = c(1,2))

with(data = subset(FV2Comp, Landrace == "KE"),
       plot(y = DH_FV2Comp, x = MinSubgroup, type = "b",
            main = "DH_KE", ylab = "FV2 Composition", xlab = "Min Subgroup"))
  
with(data = subset(FV2Comp, Landrace == "PE"),
       plot(y = DH_FV2Comp, x = MinSubgroup, type = "b",
            main = "DH_PE", ylab = "FV2 Composition", xlab = "Min Subgroup")) 

sel <- with(CV_sum, (Scenario == "Scenario2") &
              ((HBMethod == "FixedHB" & Param1 == "win20")|
                 (HBMethod == "HaploView" & Param1 == "GAB")|
                 (HBMethod == "HaploBlocker" & grepl(Train, pattern = "KE") &
                    Param1 == "MinSubgroup40" & Param2 == "MCMB1")|
                 (HBMethod == "HaploBlocker" & grepl(Train, pattern = "PE") &
                    Param1 == "MinSubgroup40" & Param2 == "MCMB1250"))
            )

dat_plot <- subset(CV_sum, 
                   (HBMethod == "HaploBlocker" &
                      grepl(Train, pattern = "KE") &
                      Param1 == "MinSubgroup30") |
                     (HBMethod == "HaploBlocker" &
                      grepl(Train, pattern = "PE") &
                      Param1 == "MinSubgroup50")
                   )
dat_plot_snp <- subset(CV_sum, 
                       Scenario == "Scenario2" & HBMethod == "SNP")

key <- match(paste(dat_plot$Trait, dat_plot$Train),
             paste(dat_plot_snp$Trait, dat_plot_snp$Train))
dat_plot$PA_SNP <- dat_plot_snp[key, "PA"]
dat_plot$HBMethod <- factor(dat_plot$HBMethod,
                              levels = c("FixedHB", "HaploView", "HaploBlocker"))
  

# Fig_scen2_SNPvsHB #  
ggplot(data = dat_plot, aes(x = PA_SNP, y = PA, color = Train, shape = Trait)) +
    geom_abline(slope = 1, intercept = 0) +
    geom_point(size = 3) +
    scale_shape_manual(values = c(8, 15:18)) +
    facet_grid(cols = vars(HBMethod)) +
    theme_bw(base_size = 20) +
    labs(title = "Scenario 2",
         x = "Prediction Accuracy (SNP)",
         y = "Prediction Accuracy (HB)",
         color = "Population", shape = "Trait") +
    scale_color_manual(values = c(wes_palette(name="Darjeeling1")[c(1,3,5)],
                                  wes_palette(name="Rushmore")[c(3)]))   

Weighted GRM

dat_plot <- subset(CV_sum, Scenario == "Weight")
dat_plot$Scale <- as.numeric(dat_plot$Scale)
dat_plot_snp <- subset(CV_sum, Scenario == "Scenario1" &
                         HBMethod == "SNP" &
                         Train == "DH_PE")
dat_plot_snp$Weight <- NULL; dat_plot_snp$Scale <- NULL


ggplot(data = dat_plot, aes(x = Scale, y = PA)) +
  geom_point() +
  geom_line() +
  geom_hline(data = dat_plot_snp, aes(yintercept = PA), linetype = "dashed") +
  facet_grid(rows = vars(Trait), cols = vars(Weight), "free") +
  theme_bw(base_size = 20) +
  labs(title = "Scenario 1 (DH_PE)",
       x = "Scaling Factor (s)",
       y = "Prediction Accuracy") +
  scale_x_continuous(breaks = seq(0, 2, by = 0.5))